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ABSTRACT 

When protoplanets growing by accretion of planetesimals have atmospheres, 
small planetesimals approaching the protoplanets lose their energy by gas drag 
from the atmospheres, which leads them to be captured within the Hill sphere 
of the protoplanets. As a result, growth rates of the protoplanets are enhanced. 
In order to study the effect of an atmosphere on planetary growth rates, we 
performed numerical integration of orbits of planetesimals for a wide range of 
orbital elements and obtained the effective accretion rates of planetesimals onto 
planets that have atmospheres. Numerical results are obtained as a function 
of planetesimals' eccentricity, inclination, planet's radius, and non-dimensional 
gas-drag parameters which can be expressed by several physical quantities such 
as the radius of planetesimals and the mass of the protoplanet. Assuming that 
the radial distribution of the gas density near the surface can be approximated 
by a power-law, we performed analytic calculation for the loss of planetesimals' 
kinetic energy due to gas drag, and confirmed agreement with numerical results. 



- 2 - 



We confirmed that the above approximation of the power-law density distribu- 
tion is reasonable for accretion rate of protoplanets with one to ten Earth-masses, 
unless the size of planetesimals is too small. We also calculated the accretion 
rates of planetesimals averaged over a Rayleigh distribution of eccentricities and 
inclinations, and derived a semi-analytical formula of accretion rates, which re- 
produces the numerical results very well. Using the obtained expression of the 
accretion rate, we examined the growth of protoplanets in nebular gas. We found 
that the effect of atmospheric gas drag can enhance the growth rate significantly, 
depending on the size of planetesimals. 

Keywords: planetary formation; planetesimals; origin, solar system 



Introduction 



Planets are thought to be formed in circum-stellar gas disks by accumulating a large 
number of planetesimals. When a protoplanet reaches about the Moon mass, it starts to 
attract surrounding gas of the disk by its gravity, which would lead to the formation of a 
thick atmosphere, if the growth of the protoplanet proceeds in the nebular gas. Thus, in 
the process of planetesimal accretion onto planets, interactions between planetesimals and 
atmospheres should be considered. Interactions between planetary atmosph eres and plan- 



etesim als have been studied mainly in the context of giant planet formation. iPodolak et al. 



( 1l988l ) studied the interaction of planetesimals with atmospheres in order to examine the ab- 
lation and dissolution of incoming planetes imals due to gas dra g by calculating traiectories 
and mass-loss rates of the planetesimals. iPoUack et al.l (119961 ) and iHubickyj et al.l (120051 ) 
followed the method of iPodolak et al.l (119881 ) to calculate energy and mass deposition on the 
atmosphere when planetesimals are passing through it as a part of their formation model of 
giant planets. Also, in relation to the origin of natural satellites, the interaction between pri- 
mordial atmospheres and incoming objects was considered, in order to study the possibility 
of capture of satellites by the gas drag, as well as the orbital evolution of captured satellites 
(IPoUack et al1ll979l : ICuk and Bu^l2004h . 



Interactions between planetesimals and atmospheres also affect the growth rate of pro- 
toplanets. In protoplanetary disks, solid protoplanets form first, and then solid protoplanets 
with a sufficiently large mass ( about ^,10 Earth masses) acquire a large arnount of gas to 



1996: Ikoma et al 



becorne gas-giant protop l anets (|Mizundll980l: Bodenheimer and Pollack 



1986; Pollack et al. 



2000l : iHubickyj et al.ll2005l ). This scenario is widely accepted, but the 



growth timesca le of sohd planets predicted by results of th eoretical studies of planetary 



accretion (e.g., iTanaka and Idal Il999l : iKokubo and Idal l2000l ) are too long in two senses 



3 



First, the predicted growth timescale is longer than the typical timescale of in ward migra- 



tion o: 



1997 



protoplanets predicted 



3V cur r ent theories includin g linear analyses (e.g.. IWardlll986 



Korvcanskv and Pollack 



D'Angelo et al 



2002 



1993 



Masset et al 



Tanaka et al.l |2002| ) and numerical simulations (e.g. 



20061 ). which results in the loss of forming planets into 
the central star before the completion of their growth. Second, planets cannot reach a critical 
mass that is necessary for the onset of runaway gas accretion before dissipation of disk gas, 
which means that gas giant planets are difficult to be formed. Curre nt theories predict that 



solid planets grow in proto-planetary disks via runaway growth ph ase (jWetherill and Stewart 



19891 ) followed by ohgarchic growth phase flKokubo and Idalll996l ). and in both phases, pro- 
toplanets grow mainly by capturing objects that are much smaller than the protoplanets. 
Thus, in order to estimate the growth timescale of solid planets, one has to investigate the 
stage where a large number of small objects are accreted by a small number of large objects. 

Collision rates of small objects onto protoplane ts have been studied in detail by numer 



ical orbital integration and analytical calculations 


Nishida 


19831; 


Wetherill and Cox 


1985 


Ida and Nakazawa 


1989; 


Greenzweig and Lissauer 


1990. 


I992I: 


Dones and Tremaine 


1993) 



The early stage of planetary accretion is likely to have proceeded in the presence of the 
protoplanetary gas disk, but most of the previous studies did not consider the effect of at- 
mosp here around protop lanet s, which likely enh a nces effective accretion rates through gas 
drag. iKary et al.l (119931 ) and iKary and Lissauerl (119951 ) studied accretion rates when plan- 
etesimals are migrating radially inward by gas drag from the nebular gas, but they did not 
consider the effect of planets' atmospheres. 



Inaba and Ikomal (120031 ) first examined the effect of atmospheric gas drag on planetary 
accretion rates in detail. They showed that the presence of atmospheres around planets can 
enhance the accretion rate greatly depending on the size of objects that are passing through 
the atmosp heres. The results were used in their s imulation of planetary a ccretion in the Jo- 
vian r egion Jinaba et al.ll2003h . Ichambersi J2006al lbl) also used the results of llnaba and Ikomal 
( 120031 ) in his semi-analytic model of oligarchic growth of planets, and showed that the en- 
hancement of the accretion rate can shorten the growth timescale so that gas giant planets 
can be formed be fore the dispersal of the gas disk, if incoming planetesimals are small (< 
100m). Although llnaba and Ikomal (120031 ) demonstrated the importance of the atmosphere 
for planetary accretion rates, most of their results were based on calculations neglecting the 
solar gravity, and three-body orbital integration was performed only in the case of initially 
circular orbits. In addition, they used a realistic model for the atmospheric structure rather 
than a simplified one, which makes the study of parameter-dependence or comparison with 
analytic calculations rather difficult. 



In the present work, we obtain accretion rates of planetesimals onto protoplanets that 
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have atmospheres, by using three-body orbital integration for a wide range of parameters. 
We use a simple power-law density distribution model for the atmosphere, which allows a 
systematic study of the accretion rate and its parameter-dependence. From the obtained 
accretion rates, we derive a semi-analytical expression of the accretion rates for planets with 
atmospheres and examine growth of protoplanets in the nebular gas. In ^ we first show the 
basic formulation and some analytic estimations which are the basis of our numerical simu- 
lation and analysis to obtain the accretion rates. In ^ we show the results of our numerical 
simulations to obtain the accretion rates, together with some demonstrative calculations, 
which are useful in understanding the basic behavior of planetesimals' motion under the 
gas drag of planetary atmospheres. In ^ we derive a semi-analytic formula for accretion 
rates under the assumption of power-law density profile for atmospheres. In ^ we examine 
the validity of the power-law approximation for density profile of atmospheres. Using the 
obtained semi-analytic formula for accretion rates, we examine the growth of protoplanets 
with atmospheres in ^ Our conclusions and discussion are presented in ^ 



When the masses of planetesimals and the planet are much smaller than the solar mass 
and their orbital eccentricities and inclinations are sufficiently small, their motions in a 
rotating coordinate system can be described by linearized equations, called Hill's equations. 
We use a coordinate system centered on the planet, where the a;-axis points radially outwards, 
the ?/-axis points in the direction of the planet's orbital motion, and the 2;-axis is normal 
to the x-y plane. We scale the distance by the mutual Hill radius rn = ah with h = 
{(M + m)/3M4i/3 (TVf and m are the masses of the planet and a planetesimal, respectively; 
a is the semi-major axis of the planet; and is the mass of the central star) and the time 
by Vt^ (^K is the planet's orbital angular frequency). The non-dimensional equation for the 
relative motion between the planet and a planetesimal can then be written as 



2. Basic Formulation and Analytical Estimate 



2.1. Basic equations 



dv 
di 




where v is velocity, t is time, and $ is the Hill potential given by 



~ 3 3 9 X ^9 9 

$ = X + -z'^ + - 

f 2 2 2 



(2) 
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where f = -^/a;^ + + 5^ . Variables with tildes denote non-dimensional quantities. The 
acceleration due to gas drag ajrag is described by 

Fdrng/m 3 Pg^_i , , 

Ctdrag = = -o'-^D^s AuAlt, (3) 

where -Fdrag = —{C-o /2)Tirlp^/S.u/S.u is the drag force for a planetesimal with radius rg, 
m is the mass of the planetesimal, Cd is the gas drag coefficient of order unity, Pg is gas 
density, ps is the internal density of planetesimals, fg is the normalized physical radius of the 
planetesimals, and Ait is the velocity of the objects relative to the gas. 



2.2. Atmospheric Structure 



In the present work, we assume that the density structure of the atmosphere is spher- 
ically symmetric. For sufficiently massive solid protoplanets (> 10 Earth masses), runaway 
gas accretion occurs and the density distribution of the inflowing gas significantly deviates 
from spherical symmetry. However, spherically symmetric atmosphere is a good approxi- 
mation for atmospheres of low-mass protopl anets because h ydrostatic (or quasi-hydrostatic) 
equilibrium is reahzed in such cases 



le.g. 



Mizuno 1980: Bodenheimer and Pollack 1986 



Pollack et al.l Il996l : llkoma et al.ll2000l ). In addition, in order to study the gas drag effect 
systematically, we assume a power-law function for the atmospheric structure as 



Pg = Ph 



(4) 



where pb is the density at the bottom of the atmosphere, or equivalently, the density at the 
surface of the solid core, f is the distance from the center of the planet, fp is the radius of 
the solid part of the planet, and the exponent a is assumed to be constant. 



For a purely radiative atmosphere with a constant opacity , the density distribution is 
given by the above power-law with a = 3 (e.g.. IStevensonI fjl982l ): see Appendix A). In more 
realistic atmospheric structures, on the other hand, the slope of the density distribution 
changes near the sublimation points of ice and sihcate dust, espe cially when opacity (or 
dust/gas ratio) is large (e.g., Ilnaba and Ikoma] l2003l : iRafikovl |2006| ). and outside of Bondi 
radius where gas is no longer bounded by the planet gravity. However, the overall structure in 
most radial regions of atmospheres can still be well approximated by a power-law distribution. 
In the present work, we adopt the above simple model for atmospheric structures in order 
to examine the effect of atmosphere systematically for a wide range of parameters and to 
better understand physics behind it. Effects of a more realistic atmospheric structure based 



on the analytic solution described in Ilnaba and Ikomal (120031 ) will be discussed in Section |5l 
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Once a power-law function is assumed for gas density, effects of gas drag can be described 
by two parameters: the exponent a, and the non-dimensional gas drag parameter C,, which 
will be defined in the next subsection. 



2.3. Non-dimensional gas drag parameter ^ 

Assuming that density profile is expressed by Eq. (jl]), Eq. ([3]) can be rewritten in the 

form 

fldrag = -^f'^AuAii, (5) 

where ^ is a non-dimensional parameter representing the strength of gas drag, and is defined 
by 

/- 3 „ Ph _Q 



P'. 



lAU/ V^ffi/ Vlkm/ Vlkgm^V \^l03kgm^ 



In the above, is the Earth mass, pcorc and are the internal density of planets and the 
Earth, respectively, and we assumed a = 3 as mentioned above. For a purely radiative atmo- 
spher e with large optical th i ckness , pb oc M^L~^f~^ is a good approximation (see Appendix 



A or ISasaki and Nakazawal (119901 ) ) , where L is the released energy at the solid surface of 
the planet per unit time due to accretion of incoming planetesimals, and is the opacity 
depletion factor, that is, = 1 when dust to gas ratio is the same as the ratio in interstellar 
clouds assuming dust optical property is the same as that of interstellar dust; decreases 
when dust itself is depleted or when dust opacity decreases (due to growth of dust, for ex- 
ample). Thus Pb, which is a function of the planet mass, can be eliminate d from the above 



equat ion. Using the analytic formula for atmospheric structures given by llnaba and Ikoma 



(120031), we found that pb can be approximated by pb ^ Q.n]<.gTnr^{M /M^f{L/l{)-'^ Lq)'^ 
(Pcoro/Pffi) (/k/0.01)~^ {Lq is the solar luminosity) when is around 0.01. Substitution of 
this Pb into Eq. ([6]) yields 



2.4 X 10"^ 




M \ f \-W p., 



lO-^Mgj/yr / Vlkm/ VlO^kgnr 



0.01 



(7) 
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where we used accretion energy L = GMM/vp and M is the accretion rate of planetesimals 
onto tlie planet. 



2.4. Normalized accretion rate 



Following the definition of collision rate given by iNakazawa et al.l (119891 ). we define the 



accretion rate of planetesimals per unit surface number density as 

^acc(e,i)= / p^ccie,i,b)^\b\db, 



where e and i are planetesimals' orbital eccentricity and inclination divided by h, and b is 
the initial value of the normalized semimajor axis of an incoming planetesimal relative to 
the planet at y = oo. Pacc(e, ^,6) is the accretion probability for a given set of the orbital 
elements {e,i,b) defined by 

~~ f f ,~~7 .drdzu 
Pacc{e,i,b) = J J 9?acc(e,«,0,r,ro)-^^--p-, (9) 

where r and w are the initial phase angles for epicyclic motion in horizontal and vertical 
directions, respectively, and y9acc(e, 6, r, tu) is a judgment function whether an object is 
accreted by the planet: unity if the object is accreted, and zero otherwise. In general, 
accretion can happen in two ways: (i) a planetesimal hits the planet regardless of whether 
it loses enough energy to become bound, or (ii) a planetesimal loses enough energy through 
gas drag to become gravitationally bound. We thus define the following two quantities in a 
similar way: Pcoi is the normalized accretion rate due to direct collision with the solid surface 



of a pl anet without an atmosphere, which is the definition identical to that of iNakazawa et al. 
( 119891 ): and Pcap is the normalized rate of capture of planetesimals within the Hill sphere of 



the planet due to gas drag, with the assumption that the planet is a point mass (i.e., direct 
collision onto the planet is neglected) . We also define pco\ and pcap in a similar way. 

When a planet does not have an atmosphere and the random velocity of planetesimals 
is large enough, the rate o f direct collision onto the sohd s urface of the planet with radius 
fp is analytically given by (jCreenzweig and Lissauer 1990 ) 



col 



^ ~ ~ 3 

Pcoiie, i,b)-\b\db 



(10) 
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where K and E are the complete elhptic integrals of the first and second kinds, respectively, 



with k = y 3e^/4(e^ + i^). Collision probability pcoi can be written as 



K / g^ + z^- (3/4)62 , 6/fp , 
Pco\[e,i,h) = — ^ 1 + — ^- -— ~- . (11) 



37r|6|iV e2-62 V + - {?> / A)h 



2.5. Analytic estimate of accretion rates due to atmospheric drag 

Without energy dissipation, objects coming into the Hill sphere always escape from the 
sphere in the end unless they collide with the planet, but they can be ca ptured within the 



Hill sp here if there is energy dissipation such as atmospheric gas drag. Ilnaba and Ikoma 



(120031 ) introduced "enhanced radius" of a planet due to atmospheric gas drag; planetesimals 
which come closer than this radius lose a significant amount of kinetic energy and cannot 
escape out of the Hill radius, thus it can be regarded as an effective radius of the planet for 
capturing planetesimals. In the present study, we assume a power-law density distribution for 
atmosphere, which enables us to obtain the expression for the enhanced radius analytically. 

We first e stimate energy dissipatio n of objects that pass through the atmosphere in a 



way similar to Ilnaba and Ikomal (120031 ). We assume that most of energy dissipation in one 
encounter occurs near the point of closest approach to the planet, because the gas density 
increases steeply with decreasing distance from the planet. We will confirm the validity of 
this assumption later, using orbital integration ( §3.2p . Energy dissipation due to gas drag 
A£' is then estimated by the work acting on a unit mass of the object: 

= /cer^rfe + z^V (12) 



where adrag = |a^drag|; is the path length of the orbit near the point of closest approach 
and is substituted by fcfmm (/c is a correction factor of order unity, and fmin is the mini- 
mum distance from the origin); and Wqo is the relative velocity between the planet and an 
approaching planetesimal before acceleration due to the planet's gravity. When the relative 
velocity is dominated by the random velocity rather than Kepler shear, can be approxi- 
mated by the relati ye velocity at the orig in in the unperturbed, non-gravitating solution to 



the Hill's equation (INakazawa et al.lll989l ) given as 
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The total specific energy of a particle on the Hill coordinate system is given by 

1 



E 



2^ '8 2' 



(14) 



wher e w is the velocity of the particle, and E corresponds to the Jacobi integral (INakazawa et al 
19891 ). Note that the orbital elements e, i, and b in the last expression of Eq. (fT^ are the 
initial unperturbed values, thus the term for mutual gravity (— 3/f) is neglected in this ex- 
pression. Since the $ = contour surface is closed, a planetesimal cannot escape from th e 
planet's Hill sphere, if the total energy E becomes negative due to gas drag (IOhtsukilll993l ). 
Therefore, the condition for capture by gas drag is given bj0 



AE > -vl + 



9 



(15) 



The critical value of the minimum distance fmin obtained by equating both sides of Eq. f|T5!) 
is denoted by Rc, which can be regarded as the enhanced radius of a planet insid e which 
incoming objects get effectively captured by the gas drag (jinaba and Ikomal l2003l ). Thus, 
using Eqs. fll2l) and f[T^ . the enhanced radius is obtained as the solution of 



vL + 9 ' 



0, 



(16) 



and is shown in Fig. [T]as a function of Voo- Depending on Voo, Rc can be divided into three 
regimes and approximated by 



(4/ce/3)^/" 
Rr ^ < (12/,^f;-2)V" 
(2/^^)i/(«-i) 



for < 9, 

for 9 < e < 6(2/ce)-'/("-'\ 
for 6(2/^0" 



;i7) 



-l/(a-l) ^ ^2 



In the low velocity regime, the effect of solar gravity is important. In the interrnediate 



Greenzweig and Lissauer 1119901 ) 



velocity regime, the two-body approximation is valid (e.g 
and Voo is smaller than the escape velocity aX r = Rc as well. In the high velocity regime 
Voo is larger than the escape velocity at f = 



i), namely AE > (l/2)v^ + 3 



^The corresponding condition given by llnaba and Ikomal (|200 
notation, is slightly different from Eq. p^ . because they did not consider the contribution of the tidal 
potential. 
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Once we obtain the enhanced radius, the capture probabihty pcap can be obtained by 
replacing fp in Eq. ffTTl) with Rc, namely 

~. 2Rl /e2 + 22 _ (3/4)52 / 6/^^ 

37r|6|z V e2-62 \^ e2 + z2 _ (3/4)52 

Then, the capture rate Pcap can be obtained by integrating Eq. (fT8l) with respect to 6, as in 
the case of Pcoi (Eq. lITUi) ). Note that the integration cannot be done analytically because 
Rc is not a simple function but a solution of Eq. f|T6|l . and is a function of h through 
(Eq.dH). 




3. Numerical simulation 

3.1. Numerical Method 

We integrate Eq. ([T]) for particles with various initial orbital elements, using the eighth- 
order Runge-Kutta integrator. For a given (e, -i), we calculate a phase volume for accretion 
in the (6, r, vj) phase space by finding orbits that lead to accretion. For high velocity cases, 
the phase volume is much smaller than the total volume, and we need to calculate a large 
number of orbits to obtain Pace with sufficient accuracy. We thus adopt a kind of adaptive 
mesh refinement method for r and w mesh; phase volume for accretion is narrowed down by 
several steps of calculations fr om a coarser mesh with a larger target to a finer mesh with a 



smaller target (jOhtsukilll993l ). In our numerical code, the initial azimuthal distance is set 
to yo = 50; we confirmed that this is large enough to obtain accretion rates with sufficient 
accuracy. We stop our orbital integration when one of the following conditions is met: (i 
f > yo + ^0 + 2eo + iq (the subscript denotes initial values), (ii) f < fp, (iii) P < 
Hereafter, we set a = 3 unless otherwise stated. 



Gas velocity is assumed to be zero on the rotating coordinate system, i.e.. An = v where 
V is the velocity vector of particles, since the gas velocity at a distance f from the planet 
is much smaller than the Keplerian velocity around the planet (= \/S/f) at the location, 
which is the typical velocity of a planetesimal moving around the planet under the strong 



^ To be captured by the protoplanet's Hill sphere, one may think that the third condition should be i5 < 
and f < 1 as well, but we consider the range of initial orbital elements corresponding to £' > because 
planetesimals cannot enter the Hill sphere if i? < initially, and also we set the gas drag term effective only 
within the Hill sphere, that is, where f < 1 and $ < 0, and E decreases only by gas drag, thus the condition 
f < 1 is automatically met when E decreases to a negative value. 
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effect of the planet's gravity. As we consider atmospheres in hydrostatic equihbrium before 
runaway gas accretion, gas velocity is much smaller than the sound speed. Also, the sound 
speed in an atmosphere well inside of the Bondi radius is always smaller than the Keplerian 
velocity. Thus, the neglect of gas velocity is consistent with the assumption of atmospheres 
in hydrostatic equilibrium. We do not consider ablation of incoming particles when they are 
passing through an atmospher e. Although ablation may be important for atmospheres of 



massive protoplanets (> lOM^: iBenvenuto and Bruninill2008l ). its effect ca n be neglected for 



atmo spheres before runaway gas accretion considered in the present work (jinaba and Ikoma 
2003h . 



3.2. Effects of gas drag 

Before discussing numerical results of accretion rates in detail, first we show some ex- 
amples of orbital calculations, which help us understand the basic dynamical behavior of 
planetesimals affected by gas drag in atmospheres. 

Figure [2] shows examples of coplanar orbits under gas drag of atmospheres, with three 
different values of the gas drag parameter ^. In the case with ^ = 1.66 x 10~^ (left panel), 
the orbit is slightly affected by gas drag, but the energy dissipation was not significant and 
the object escapes from the Hill sphere. When = 1.67 x 10~^ (middle panel), the particle 
loses enough energy by gas drag to be captured within the Hill sphere. When ^ = 2 x 10~^, 
the particle loses a large amount of energy at the first encounter, and spirals toward the 
planet quickly. 

Next, we examine the change of energy of a particle. Figure [3] shows the change of E for 
the captured orbit described in the middle panel of Fig. [2l The left and right panels show 
the plots of as a function of time and the distance from the planet, respectively. We find 
from the right panel that the energy reduces greatly when the particle passes through the 
dense part of the atmosphere near the planet's surface (i.e., f < 0.1 in the case shown in Fig. 
[3]), whereas it is almost constant when it is far from the planet. Also, the left panel shows 
that significant energy decrease occurs instantaneously, because most of energy dissipation 
occurs near pericenters, where particles move with the fastest velocity in an orbit. Note 
that, for the purpose of demonstration, the condition (iii) for stopping integration described 
above was not applied to the orbital integrations shown in Figs. [2] and [31 

Figure H] shows the total change of energy E due to gas drag through an encounter with 
the planet until one of the truncation conditions for orbital integration (i)-(iii) is met, as a 
function of the minimum approach distance from the planet's center. Two groups of points 
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with AE > 1 X 10~^ (red pluses and green crosses) show numerical results for the low- and 
high-velocity cases, respectively. We can see that the energy dissipation is well described by 
the analytic estimates given by Eq. f[T^ . which are shown by the two lines (we set = 2 
in Eq. flT^ ). For the low velocity case, AE oc r~^^, which corresponds to the second term of 
Eq. (fT2|) and means that the particles are well accelerated by the planet gravity at f ~ rmin- 
For the high velocity case, AE oc rl^^, which corresponds to the first term of Eq. (fT2|) and 
means that the particle velocity around f ~ fmin is not much larger than Voo- We also plot 
the changes of E for calculations without gas drag, which corresponds to the numerical error 
in our numerical integration. We find that the numerical error is much smaller than the 
calculated energy changes due to gas drag, even in the case when gas drag is very weak 
(e = 1 X 10-9). 

In Fig. [5l we show the total change of E in the same way as Fig. H] in the case with 
fp = 0.005 to examine the dependence on ^ and e. Note that there are cases when the second 
or a later close approach reduces energy by an amount comparable to or larger than that 
due to the first one even when planetesimals are not captured within the Hill sphere yet, but 
such events are rare unless planetesimals are captured by the first close approach. Thus, the 
values of AE shown in Fig. [S] practically represent the energy change due to the first close 
approach to the planet, which can be used as an indicator whether planetesimals are captured 
or not from the point of view of energetics. The left panel shows the dependence on ^ for a 
low-velocity case. In this cases, we find AE oc f"^ for most points, as explained above. In 
the weak gas drag case, planetesimals reduce their energy due to single close encounter and 
escape from the Hill sphere. In the strong gas drag case, points shift upward in accordance 
with stronger energy dissipation. The points in this case fall on a horizontal band with 
AE ~ 2-3 at fmin < 0.015; in this regime, particles are captured by gas drag rather than 
direct collision onto the planet, i.e., orbital integration is stopped by the condition E < 0. 
Thus, values of AE in this horizontal band in the plot show the amount of energy that 
needs to be dissipated for capture, and the upper end of this horizontal band (rmin ~ 0.015) 
corresponds to -Rc (see Eq. f|T7|) ). When gas drag is weaker, energy dissipation becomes 
smaller and Rc becomes smaller accordingly, and eventually the trend AE oc f-^^ does not 
reach the energy needed to be captured [AE ~ 2-3 in this case) before it reaches the surface 
of the planet (f = fp); in this case, particles are accreted by direct collision. Note that 
the points with fmin ^ 0.005 represent cases where orbital integrations were stopped by the 
condition of direct collision to the planet. The right panel shows the dependence on random 
velocity (i.e., e and i). The low velocity case is identical to the strong drag case in the left 
panel. With increasing random velocity, the energy needed to be captured increases because 
initial kinetic energy is larger, and Rc decreases accordingly. In the high velocity case, the 
numerical results fall on a line with a nearly constant slope which extends all the way to 
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'"min = fp, which Hieans that gas drag is not strong for particles with such a high random 
velocity and they are accreted only by direct collision with the planet's solid surface. 



3.3. Accretion rates for the case with e = 2i 

Next, we examine accretion rates in the case of e = 2i. Figure E] shows the plots of 
Pcoi{e,i,b) and Pcap(e, as a function of b, for two different sets of (e, z). Crosses and 
pluses represent numerical results for pcoi (Eq. ffTT]) ) and pcap (Eq. ffTSl) ). respectively, and the 
dotted and solid lines show the analytic results. Note that we empirically set fc = 1.83 for 
the analytic results of Pcap! changing only shifts the lines in the vertical direction and does 
not change the functional form. In the high (e, i) case, we find that the numerical results are 
in good agreement with the analytic solutions over a wide range of b, except for the small-fe 
region where some small deviations are observed. In the low (e, i) case, numerical results 
deviate from the analytic solution significantly. This is because Wqo is not large enough to 
neglect the three-body effect in this case. 

Figure [7] shows numerical results of Pace for several values of ^, along with Pcoi (bottom 
line). We also plot numerical results for PcoLhui, which is the collision rate onto the Hill 
sphere (top line). All these curves can be roughly divided into two parts: a horizontal part 
in the low-velocity regime and a part with a negative slope in the high-velocity regime. In 
the lower eccentri city regime, previous studies showed that P mi is independent of Voo and 



depends on Vp as (jida and Nakazawa Ill989l : llnaba et al.ll200ll ) 



col, low 



ll-S^rp. (19) 



This corresponds to the horizontal part of the curve for Pcoi in Fig. [3, even though the 
horizontal part can be seen only in the lowest e regime displayed in the left figure. In the 
high eccentricity regime, Pcoi can be well described by Eq. ffTOl) and is thus proportional to 
e~^; this corresponds to the negative-slope part of the above curve. 

In the case for planets with atmospheres, this general trend still holds, but changes in 
some parts. In the lower eccentricity regime, Pcap is expected to be proportional to ^^/'^"'^ 
because Pcap oc R^J'^ (see Eq. ^) and ~ (4/c^/3)^/" (see Eq. On the other 

hand. Fig. [7] shows a rather weak dependence on ^ in the low-velocity regime, especially 
when ^ > 10~^. In order to examine the ^-dependence in the low-velocity regime in detail, 
we calculated capture rates as a function of ^ for e = i = (Fig. [8]). In the low ^ regime 
~ when a = 3), Pcap is indeed oc ^^Z^". in the high ^ regime > 1), on the other 
hand, Pcap is constant and is equal to the collision rate onto the Hill sphere, because particles 
entering the Hill sphere immediately get captured due to the strong gas drag. In between 



-14- 



the above two extreme regimes of ^ in Fig. [HI there is a part where Pcap is well approximated 
by at 10~^ ^ ^ 1 when a = 3. We thus obtain the following empirical formula for 

capture rates in the low velocity regime 

-fcapO = niin(-PcapO,low,t; -PcapO.mcdg; -fcolO,Hill) ; (20) 

where Pcapo,iow$ is the capture rate in the low ^ regime, -Pcapo.med^ is the one for the interme- 
diate ^ regime, and Pcoio.Hiu is the collision rate onto the Hill sphere, respectively, and these 
variables are defined in the case with e = i = 0, given by 

-PcapClow? = 10.7 ( ^^-^ j , (21) 

-PcapO,modC = -fcol.HillC"^^^^" , (22) 

Pcoio.Hiu = 4.392. (23) 

As we can see in Fig. [8|, the expressions fl20l) - fl23l) reproduce our numerical results very well, 
even when a = 2 and 4. 

In the high-velocity regime, we can see from Fig. [7|that, when e ~ 10 and is small, 
the e-dependence of Pcap is stronger than that of Pcoi (oc e~^), because Rc decreases with 
increasing random velocity in this regime (see Eq. (ITTl) ). Thus, with increasing e. Pace for 
small ^ approaches Pcoi at a certain e (e ~ 5 when ^ = 10"'' and fp = 0.005, for example). In 
^ we will derive an analytic expression of the capture rate in this regime. The expression 
shows that the capture rate is expected to be proportional to e~^/^ when a = 3 and e = 2i, 
which agrees well with the above numerical results. 



3.4. Accretion rates averaged over the distribution of eccentricities and 

inclinations 

Planetesimals in protoplanetary disks have a distribution of e ccentricities and inclina - 



tions that is well described by the Rayleigh distribution function (jida and Makino Ill992l ). 
thus we next show Pace for incoming planetesimals that have the Rayleigh distribution in 
eccentricity and inclination. 

Normalized accretion rate for planetesimals that have a Rayleigh distribution in e and 
i with the root mean square values e* and i* is given by 

(Pace) = [[ Paee/R(g,^,g^^*)«, (24) 
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where /r is the Rayleigh distribution function given by 

4ez 



/R(e,i,e*,r) 



exp 



e2 



;*\2 



(25) 



In order to obtain (Pace) from Eq. fl24|) for certain regimes of e* and z* (10"^ < e* < 10^/^ 
and 10~^ < < 10^/^ in this paper), we need value s of P^f^n for a wide range of e and i 
including the region outs ide of the regim e of e* and i* flOhtsukilll999l ). To do that, we follow 
the method described in lOhtsukil (119991 ) with a slight modification. 



We first obtain Pacc(e, i) by direct orbital integration described in the previous sections 
at the 10 X 10 mesh points in the square region 10~^ < e < 10^^^ and 10"^ < < 10^/^ with 
even interval in the logarithmic space. We also perform orbital integration at supplementary 
coarse grid points in the low inclination regime, i.e., 3x2 mesh points with e = 10"^, 10", 10^ 
and i = 10~^, 10~^. Values outside of the region are obtained by extrapolation from the 
calculated region and by the analytic estimation for the high- velocity regime (see We 
confirmed that (Pace) for 10~^ < e* < 10^/^ does not depend sensitively on the Pace values 
for the outside of the square region 10-1 <e< 10^/^ and 10"^ < i < 10^/^ 

Symbols in Fig. [9] shows numerical results for (Pace) as a function of e* when e* = 2i*. 
While this figure looks similar to Fig. [TJ the results shown in Fig. [9] change more smoothly 
with increasing random velocity than those of Fig. [71 and the absolute values in Fig. [9] tend 
to be slightly higher than those of Fig. [7] in the high- velocity regime. These two differences 
come from the average over the Rayleigh distributi on, and similar results were obta ined for 
the Rayleigh distribution average of collision rates flGreenzweig and Lissauer Ill992l ). 



4. A semi-analytic expression for accretion rates 

Next, we derive a semi-analytical expression for (Pace) based on the numerical and 
analytical calculations for accretion rates described above. As shown in Figs. [7] and [HI 
with decreasing strength of gas drag (i.e., accretion rates decrease and approach the 
collision rates onto the solid surface of the planet without an atmosphere. Thus (Pace) can 
be approximated by the larger of (Pcoi) and (Pcap): B 

(Pace)(rp, e, «, e\ n = max((Peoi)(fp, r , r), (Pcap)(e, «, e\ n)- (26) 



■^Note that, since the integrands in the definition of Pcoi and Pcap (see Eq.®) depend on b, the 
magnitude relation between the two also depends on b, so a more precise procedure would be to take 
Pace = max (pcoh Pcap) fifst and then integrate it with respect to b, e, and i. However, for the derivation of 
an approximate formula for (Pace), the simple procedure described here is sufficient. 
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Collisi on rate averaged over the Raylei gh dist ribution was studie d in detail by previous 
work (e.g., iGreenzweig and Lissauer Ill992l ). and llnaba et al.l (120011 ) derived the following 
semi-analytic expression for (Pcoi): 

(Pcol) = min nPcol)mcd, {(Pcol)low + (-Pcol)high} ^ 



where 



( Pcol ) low 
(Pcol) med 



11.3A/f 



p 



Airi* 



17.3 



232 



(Pcol) 



high 



p 

2tv 



(27) 

(28) 
(29) 

(30) 



In the above, Eq.(l30l) was obta ined by averaging Eq. ( fTOj) over the Rayleigh distribution 
(IGreenzweig and Lissauer Ill992l ): I* = i*/e* and the function JF(J*) and Q{I*) are given by 



j^(r) 







'[(P)2+(1-(P)2)A212^ 



g{r) = 8 dx 



AV3(1 - A2)/2] 



When /* 



'[(P)2 + (1-(P)2)A2]' 

1/2, J^(l/2) = 17.34 and ^(1/2) = 38.22. 



(31) 
(32) 



As for the capture rate averaged over the Rayleigh distribution, we found that the 
following expression reproduces our numerical results very well: 

(Pcap) = ((-Pcap)iow + (-^cap)higli) ' (^3) 

with 7] = 2/3. The accretion rate for the low-velocity regime can be simply written in terms 
of Peapo i^q.m) as 



(Pcap) 



cap / low 



P 



capO ; 



(34) 

because Pcapo is independent of e and i. We further divide (Pcap)high into two functions as 

(Pcap)liigh — (Pcap)highl ~l" (Pcap)higli2- (35) 

We define that (Pcap)highi is for a relatively low-velocity regime in the sense that the relative 
velocity is low enough for the collision cross-section to be enhanced by gravitational focusing, 
whereas it is still high enough to hold the two-body approximation. On the other hand, 
(Pcap)high2 is defined for a higher velocity regime where gravitational focusing is negligible. 
In the former regime, Pcap,highi can be obtained by the integration of Eq. (ITSl) as in Eq. (jH]): 

12(12/,0'/" 



P 



cap,highl 



7rz(e2 + z2)7 



■P(A;,7), 



(36) 
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where 7 = (a + 2)/2a, and R{k,'~f) is an integral similar to the complete elliptic integral, 
and is defined by 

"^^^'^^"io (1-x^mi-Px^)^ - ^'^^ 

The Rayleigh distribution average of Pcap,highi can be calculated as 

where F is the gamma function, 2-^1 is the hypergeometric function, and we assume e* > i* 
and a > 1. This gives (Pcap)highi oc (g*)-(27+i) = (^g*)-2{a+i)/a ^j^^^^ ~* ^ 2!*, whose 
dependence on e* (i.e., oc (e*)~2("+^)/") agrees with that of Pcap.highi on e when e = 2i. For 
the typical case of a = 3 and e* = 2^*, we have (Pcap)highi — 1-52 x 10^(/cO^^^(e*)~^/^. 
In deriving Eq. (!38|) . we treated R{k,'y) as a constancl, because i?(/c,7) is a slowly varying 
function of k as long as is not close to unity (note that < A; < ^/S/A from its definition). 
For the still higher velocity regime, 

Pcap,high2 = -^-^ J{e/zy + lE{k), (39) 

thus 

(i^cap)high2 = ^-^^ HH- (40) 

The curves in Fig. [9] show (Pacc)(?^p, ^, a, e*, i*) obtained using the above semi-analytic 
expressions. We can see that (Pace) is well approximated by the above semi- analytical ex- 
pressions within a factor of two for the entire regions of random velocity and ^. When fp is 
large, (Pace) tends to approach (Pcoi) because planetesimals tend to collide with the planet 
more easily before losing a sufficient amount of energy by gas drag to get captured. For 
example, when fp = 0.005 (left panel) and ^ = 10~^, (Paee) = (Peoi), which can be deduced 
by Fig. [H On the other hand, (Paee) can be approximated by (Peap) when ^ is large. In 
summary, the effect of atmosphere is statistically negligible in the entire velocity regime 
when ^ is so small that (Pcoi)high > (Peap)highi is satisfied, where (Pcoi)high and (Pcap)highi are 
given by Eqs. fl30|) and fl38|) . respectively. This condition can be written in terms of C, as 

e < 4 X 10"^ ( . (41) 
To derive the condition, we set e* = 2 as the lowest limit of the range for the functions. 



^We confirmed that the error arisen from tliis simplification is about 3% when e* = 2i* and a = 3. 
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It should be noted that when the mass of planetesimals are small (< lO^'^g), the as- 
sumption of e* = 2i* breaks down if protoplanetary gas disk still exists a nd runaway bodies 
significantly perturb the orbit of such small bodies (lOhtsuki et al.ll2002l ). This is because, 
once inclination is damped to z < 1, the stirring rate of inclination by the large bodies is 
much smaller th an that of eccentricity in such a velocity regime, and only eccentricity tends 
to be enhanced (jldaj[l990t). 

Using analytic calculation of planetesimals' energy d issipation due to g as dr ag under 
the two-body approximation neglecting the solar gravity, llnaba and Ikomal (120031 ) derived 
a relation between the size of planetesimals and the enhanced radius of protoplanets, for 
a given atmospheric structure (their Eq.(17)). Substituting thi s enhanced radius for the 
planetary radius in the expressions for collision rates derived by llnaba et al.l (120011 ). semi- 
analytic expressions of accretion rates for protoplanets with atmospheres can be obtained 
( jinaba and Ikomal l2003l . their Eqs.(20)-(24)). We find that (Pam) with the e f fect o f atmo- 
spheric gas drag obtained from these expressions derived by llnaba and Ikomal (120031 ) ((-Pcoi) 
in their notation) is smaller by a factor of about two as compared to ours in the high- velocity 
regime, even when the same atmospheric structure is assumed (Fig JTOl) . This discrepancy 
i s cau sed by the more simplified procedures in deriving the expression in llnaba and Ikoma 
( 120031 ) as compared to ours; th e cause for the discrepan cy is further discussed in detail in Ap- 
pendix B. On the other hand, llnaba and Ikomal (120031 ) obtained the enhanced radius using 
more realistic atmospheric structures, while we assumed a power-law density distribution. 
Effects of this assumption will be discussed in Section 



5. Comparison with the case of more realistic atmospheric structures 

Here, we examine the validity of our assumption of the simplified atmospheric structure, 
i.e., power-law function (Eq.(|l])). As shown above, the assumption of power-law approxi- 
mation for gas density in atmospheres allowed a systematic study and understanding of the 
dependence of the accretion rate on various physical parameters. In order to check the va- 
lidity of this assumption, we compared the results shown previously in this paper with those 
ob tained using a more reali stic atmospheric structure based on the analytic solution shown 



by llnaba and Ikomal (120031 ). Although the index a should be properly chosen depending 
on the parameter values for a given atmospheric structure, we found that the use of power 
law functions is a good approximation in most cases to describe realistic atmospheres when 
opacity depletion factor is much smaller than unity (/^ < 10~^), because optically thin 
atmospheres tend to be radiative rather than convective. The assumption of small would 
be reasonable, because a large fraction of dust in protoplanetary disks in this stage would 
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be depleted after the accretion of large bodies such as planetesimals or protoplanets. In 
addition, small also corresponds to the situation when dust size is much larger than sub- 
micron size, and this would be the case for large protoplanets (M > M^), because dust in 
the atmos pheres grows quickly to reduce opacity by a factor of 100 from that of sub-micron 



size dust ( iMovshovitz and Podolakll2008l ). 



However, the power-law distribution would not be a good approximation at outer re- 
gion of the atmospheres in some cases. The density profile fitted by a power-law function 
at the bottom of an atmosphere tends to underestimate the gas density at the outer region 
when planet mass is small (M ~ O.IM0) and overestimate it when planet mass is large 
(M ~ lOM®) for typical disk parameters. In the case of small planets, however, the at- 
mosphere is not very effective in capturing planetesimals because it is not massive enough. 
On the other hand, the effect of atmospheric gas drag on accretion rates is most important 
for such large bodies a s M ~ lOM^, because they have thick and massive atmospheres 



( llnaba and Ikomall2003l ). Also, the formation timescale of protoplanets is determined by the 
later stage of accretion (i.e., when M ~ I-IOM® for ~ 5AU) where their growth timescale 
is long. Thus here we concentrate on the cases of large protoplanets, and compare accretion 
rates for atmospheres with power-law density profiles and those for more realistic atmo- 
spheric structures. As examples, we consider cases of IM^ and lOM^ planets at 5AU when 
disk gas temperature is 125K and density is 1.5 x 10~^kgm~^. We set the density of solid 
part of the planets to be 3.4 x lO^kgm"^ (fp = 0.001), which rep resents the density o f lOM® 



planets that consist of ice and rock for roughly even mass (e.g., iFortney et al.l 120071 ). 



Figure [TT] shows the density profiles derived by the method given by llnaba and Ikoma 



( 120031 ) together with fitted power-law profiles. As we can see, the fitted power-law functions 
overestimate the gas density in the outer region (f > 0.01), by a factor of up to ~ 10^ for the 
IOM0 case. As for the IM0 case, the line for the power-law distribution crosses the realistic 
one and underestimates the gas density at the outermost region (f > 0.2) because the region 



i s alre ady outside of the outer boundary of the analytic solution given by llnaba and Ikoma 



( 120031 ) and the temperature and density are constant. The density profile at the outer region 
also depends on the outer boundary condition, which will be briefly discussed later in this 
section. 

Figure [12] shows AE obtained by orbital integration as a function of fmin for the two 
density profiles for the lOM® planet for various planetesimal sizes (10 - 10 ''m) with e = 0.1 
and 10 in order to examine the dependence of Rc on rg (or ^) and e when e = 2i (the values 
of fmin for the upper end of the horizontal band for each curve in this plot corresponds to Rc, 
as we discussed in §3.2p . In the case of e = 0.1 (left panel), we can see from the shapes of 
the curves that the difference of AE directly reflects the assumed density profiles (Fig. [TTj) . 
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When planetesimals are large (vg = lO^km), Rc of the two atmospheres are almost the 
same because the two curves eventually agree where the curves reach AE ~ 2 — 3, which is 
necessary to be captured within the Hill sphere when e ~ 0.1. But as planetesimals become 
smaller, the curves of AE, which is dissipation energy per unit mass, shift upward, which 
yields two different Rc, the values for the realistic atmospheric model being smaller. When 
Ts =lkm, the difference in the value of Rc for the two atmospheric profiles is about a factor 
of three. In the case of e = 10 (right panel), energy dissipation required for capture is larger. 
Thus, in order to be captured, planetesimals have to pass through the inner dense region of 
the atmosphere, where the difference of the two density profiles are smaller, yielding better 
agreement in Rc for > 1km. 

Figure [13] shows the plots of Pace in the case of e = 2z for the two atmospheres for 
various planetesimal sizes. In the case of M = lOM^ (left panel). Pace for the fitted power- 
law density profile overestimates the accretion rate for small planetesimal sizes, as easily 
deduced from the results of Figs. [11] or [121 whereas the two values of Pace for the two 
atmospheric models are almost the same when planetesimals are large. Smaller planetesimal 
cases tend to give a larger difference in Pace because planetesimals are captured at outer 
atmospheres, where the difference of the two density profiles is larger. The values of Pace 
differ by a factor of 2-3 for 1km planetesimals, and the difference is about a factor of 15 for 
10m planetesimals. As for the dependence on eccentricity, there are two effects that shift 
Pace in opposite directions. One is the e-dependence of the amount of energy dissipation 
required for capture, which we already discussed above and the difference in Pace for the two 
atmospheric profiles becomes smaller for planetesimals with higher eccentricities (Fig. W2\i . 
This tendency can be seen in the case of rg = lO^km, where the difference between the 
results for the two atmospheric profiles becomes smaller with increasing e. The other is the 
effect of gravitational focusing by the planet. When the initial random velocity of incoming 
planetesimals is small (e ^ 1), velocities of planetesimals at Rc (-C 1) are determined by 
the planet's gravity, and thus the effective collision cross section is proportional to Rc {cf. 
Eq. f[T8l) ). Therefor^, the difference in Pace is roughly proportional to the difference of 
Rc. However, when the initial random velocity is large, planetesimals are not significantly 
accelerated by the planet's gravity at locations with r > Thus, when Rc > Pace 

tends to be proportional to the geometrical cross section, i.e., oc R^, which amplifies the 
difference in Pace- This tendency can be seen in the case of rg = 10m, where the difference 



^Note that Pace is a quantity integrated with respect to 6, which means that Pace is not directly propor- 
tional to Rc- But the dependence of pcoi or pcap on b is not large, as shown in Fig. [6l 

^Note that such a high velocity, which results in the large difference in Pace for planetesimals with = 10m 
size (by a factor of 15), is unlikely to be realized. Because, to enhance the relative velocity such that Paee is 



- 21 - 



between the results for the two atmospheric profiles becomes larger with increasing e. In 
the case of M = IM^, the difference between the values of Pace for the two atmospheric 
profiles can be explained basically in a similar way. In both IM^ and lOM^ cases, the 
difference between the results for the atmospheric profiles is rather small in the low-velocity 
regime with e < 1. One difference from the M = lOM® case is that Pace with the power-law 
model for the smallest planetesimal case (rg = 0.1m) underestimates the accretion rate when 
e < 10, which simply reflects the underestimation of the density at the outer region (see 
Fig. [TTl) . However, it should be noted that the assumption of our numerical method (i.e., 
planetesimals move on a Keplerian orbit while weakly perturbed by gas drag) breaks down 
for such small planetesimals. Their motion is expected to be strongly coupled with the gas 
flow, and the accretion rate of such small planetesimals is determined by the flow pattern of 
the nebular gas in the vicinity of a planet, rather than the structure of the atmosphere near 
the planet's surface. Therefore, the difference in the accretion rates for the two atmospheric 
proflles is practically negligible in the M = IM^ case, if we take into account the fact that 
sufficiently small planetesimals are coupled with the gas drag flow. 

In summary, the power-law atmospheric density proflle is a reasonable approximation 
for accretion rates of large (~ lOM^) protoplanets when planetesimals' size is not too small. 
In the case of IM^, the difference of the two atmospheric proflles practically does not cause 
a signiflcant difference in Pace, if we take account of the fact that the motion of planetesimals 
with very small sizes (< 0.1m) are strongly coupled with the nebular gas. 

We further note that the outer boundary condition of atmospheres is important for the 
structure of the outer region of atmospheres, where the difference between realistic density 
proflles and fltted power-law proflles can be signiflcant, although the inner dense region, 
which account for most of the atmospheric mass, is insensitive to the boundary condition. 
For example, when the temperature of the outer boundary is decreased to 50K, which is a 
typical disk temperature at 5AU for a disk that is optically thin in the vertical direction 
and optically thick in the horizontal direction, the agreement between the realistic density 
proflle and the fltted power- law function is signiflcantly improved (see Fig. [11]); in this case, 
the difference in accretion rates for the two atmospheric proflles is expected to be reduced 
greatly. Conversely, when the density at the outer boundary decreases due to the dispersal 
of the nebular gas, for example, the difference between the two proflle increases (Fig. [TT]) . 



roughly oc i?^, planetesimals need a close encounter that minimum distance to the protoplanet have to be 
closer to the enhanced radius, which means that planetesimals would be captured by the encounter. 
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6. Growth rates of protoplanets 



In order to examine the consequence of the obtained formula of accretion rates, we 
calculate the time evolution of planet mass using a simple semi-analytic model with the ac - 
cretion rate obtained above. We basically follow the calculation model of lChamberd (l2006al ). 
and also compare our results with his results. 

We consider a protoplanet embedded in a disk of single-size planetesimals, whose random 
velocity is determined by the equilibrium between damping by gas drag and excitation by 
the protoplanet. The growth rate of a protoplanet is given by 



dM 



27rSpr| 
T 



(-Pace) ; 



(42) 



where Sp is the surface density of planetesimals and T is the orbital period of the protoplanet. 
To integrate this equation, we need formulas for the surface density and (Pace)- The surface 
density is given by 

M 

Ep = Etot (43) 



27ra6rH 

where Etot is the initial surface density of planetesimals. This expression was derived under 
the assumption that the protoplanet does not migrate radially. In this case, the final mass of 
the planet becomes the isolation mass, (27ra^6fEtot)^^^/ (3M*)^/^, where 6f is the width of the 
feeding zone of the planet scaled by the Hill radius. As for the normalized accretion rate, we 
use our se r ni-anal ytic expression for the high- velocitv regime fEa.([38|)) in acc ordance with 
Chambers! (j2006al ) , who used the expression based on llnaba and Ikomal (120031 ) for the same 
velocity regime. Equation (138!) includes the gas drag parameter ^, which is a function of the 
accretion rate (see Eq.(I71)). Thus we solve simultaneous equations with respect to dM/dt 
(Eq.dlSl) with Eq.(l3H])) and ^ self-consistently. We set E^t = 100 kg m'^ and a = 5AU, 



and also assume that mean molecular weight is 2.8, material density is 2.0 x lO^kgm ^, gas 



density of the disk is 3.9 x 10 ^kgm ^, opacity is 1.0 x 10 ^m^kg ^, and a 
atmospheric gas density distribution. 



3 for the 



Figure [m^a) sh ows the growth of the protoplanet. We first calculated the growth with 
the formula used in 1 Chambers! (l2006a! ) to check our calculation, and confirmed that their 
results in the case of planetesir nals with 10km dia meter can be reproduced (thin dashed line 
in Fig. ITW a): see also F i g.8 of [Chambers! (l2006al )). As demonstrated by the semi-analytic 
calculation of 1 Chambers! (l2006al ). acceleration of the growth of protoplanets by the effect of 
atmospheric gas drag becomes notable when M > M^, because at this stage the atmospheres 
of the protoplanets become dense enough to capture 10-km-diameter planetesimals by gas 
drag. 
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On the other hand, the three thick hnes in Fig. [TWa) show results obtained using our 
new formula for accretion rates with the effect of atmospheric gas drag. We also plot the 
change of values of ^ in the course of planetary growth in these three cases in Fig. ITW b). We 
can see that the growth based on our formula of the acc retion rate for planetesimals with 
10km is faster than that obtained by I Chambers! (l2006ar) by a factor of abou t two, which 



reflects the difference between our formula and that of llnaba and Ikomal (120031 ) . as discussed 
in ^and Appendix B. In addition to the case of 10-km-diameter planetesimals, we also per- 
formed calculations with smaller planetesimal sizes (1km and 100m). The acceleration of the 
protoplanet's growth is much more significant in the case of small planetesimals for two rea- 
sons. First, in the early stage where M < 10~'^M(^, the growth rate is enhanced for the small 
planetesimal cases because their equilibrium random velocity is smaller due to strong drag 
by the protoplanetary disk gas than that in the case of the larger planetesimal case. Second, 
in the late stage where M > O.IM^, the enhancement of accretion rates by atmospheric gas 
drag is more significant in the small planetesimal case, because corresponding values of ^ are 
larger for small planetesimals as shown in Fig. [T4l(b). As a result, the growth is accelerated 
by more than a factor of 10 in the case of 100-m-diameter planetesimals as compared to 
the case of 10-km-diameter. We do not intend to draw many conclusions from our growth 
calculation presented in Fig. [TU because our model is rather simple and does not include 
other important effects such as fragmentation or migration. However, our calculation shows 
the importance of atmospheric g a s dra g in planetary acc r etion , as first demonstrated and 
emphasized by llnaba and Ikomal (120031 ) and llnaba et al.l (120031 ) , and it also demonstrates 
the usefulness of the new formula for accretion rates derived in the present work. A more 
detailed growth model, which includes fragmentation of planetesimals, non-equilibrium ec- 
centricity, a nd radial rnigration as w ell as planetesimal capture by atmospheric gas drag can 
be found in Ichambersi J2006al . boosh . 



7. Conclusions and Discussion 

In the present work, we examined accretion rates of planetesimals onto protoplanets 
that have atmospheres by analytic calculation and numerical orbital integration with gas 
drag. Assuming that the radial distribution of the atmospheric gas density can be approxi- 
mated by a power-law and that most of energy dissipation of a planetesimal passing through 
a protoplanet's atmosphere occurs near the point of closest approach to the protoplanet, 
we analytically estimated dissipation of kinetic energy of the planetesimal, and confirmed 
agreement with orbital integration. We performed three-body calculation for a large number 
of orbits and obtained accretion rates for a wide range of parameters. We also examined 
the case when planetesimals have the Rayleigh distribution in orbital eccentricities and in- 
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clinations, and derived a semi-analytic expression of the accretion rate with the effect of 
atmospheric gas drag. 

We found that our semi-analytic formula given by Eq. fl2B]) can reproduce our numerical 
results for the accretion rate of planetesimals with the Rayleigh distribution in eccentricities 
and inclinations by planets with atmospheres very well. The formula is basically described 
as a function of five non-dimensional parameters: the normalized radius of solid surface of 
planets fp, the non-dimensional gas-drag coefficient ^, the exponent of power-law function 
for the gas density profile a, and r.m.s. eccentricity e* and inclination i* of planetesimals. 
Many other physical parameters with real dimension such as semi-major axis, planet mass, 
radius of incoming particles, etc. are all reduced to ^ as in Eq.([7]), thus our results can be 
applied to a wide range of situations through the above non-dimensional parameters. 

We also performed orbital calculation with gas drag from an atmosphere with a realistic 
density distribution, and compared with the above results for an atmosphere with a power- 
law density distribution. We found that the results using these two different atmospheric 
profiles agree well with each other, for example, when the protoplanet's mass is lOM® and 
planetesimals are not too small {r^ > 1km), while the results using the power-law density 
distribution tends to overestimate the accretion rate when planetesimals' size is smaller 
and their random velocity is large (e ^ 1). We also found that the degree of deviation 
depends on the outer boundary condition of the atmosphere because the structure of the 
outer atmosphere is sensitive to the boundary. 

Using a simple semi-analytic model and the above results of accretion rate with atmo- 
spheric gas drag, we performed calculation for the growth of protoplanets. We confirmed the 
results of previous st udies that atmospheric gas drag can sign i ficantly enhance t he growth 



rate of protoplanets (jinaba and Ikomal l2003l : llnaba et al.ll2003l : IChambersI l2006al ) . We also 



found that the acceleration of the growth is significantly enhanced for small planetesimals. 
Since our semi-analytic formula is expressed in a general form in terms of non-dimensional 
parameters, it can be used in the study of planetary accretion in various situations. 
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Appendix A: Gas density distribution of a wholly radiative atmosphere near 

the planetary surface 

Equations of hydrostatic equilibrium and diffusion approximation of radiative transfer are 
given by 

— = (44) 
dr ^' 

IGo-T^ dT _ _ L ^^^^ 

3kp dr Airr'^ ' 

where P is pressure, G is the gravitational constant, M is planet's mass, r is distance from 
the planet center, p is gas density, a is the Stefan-Boltzmann constant, T is gas temperature, 
L is luminosity, and k is opacity. Eliminating r from equations (jS]) and (H5ll . we have 

dP GAnaGM. 

^ = ; — T\ 46) 

dT 3kL ^ ' 

Integrating this equation, we have 

P.}^r, (47) 

where temperature and pressure at the outer boundary (photosphere) is assumed to be much 
smaller than those of atmospheres near the planetary surface. Using the equation of state 
P = pR'T {R' is the gas constant), 

167raGM . , , 



Substitution of this into Eq. ( 1451) yields 

dT GM 



dr AR'r^ 

Integrating this, we have 



(49) 



GM , , 

4i?'r ^ ' 
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Substituting this into Eq. (H8l) . we have 



^ = T2^w- ^''^ 

Since M = (4/3)7rrpPcore5 where and pcore are radius and mean density of the planet, 
respectively, the density at the surface of the planet can be written in the form: 

n'^aPcorcG^M^ Til , . 



Appendix B: Comparison with llnaba and Ikomal (120031 ) 



As shown in Fig. [10, accretion rates derived in the pr esent work are sys t emat ically larger 
than those obtained by using the formulas derived by llnaba and Ikomal (120031 ). Here, we 
compare our results given by Eq. (l26l) with their results, which are obtained mainly by a 
simpler way (i.e., two-body orbital integration and rough estimation of dissipation energy of 
gas drag), and we examine the cause of the difference. 

Calculating the amount of planetesimals' kinetic e nergy dissipated by gas d rag under the 
two-body approximation neglecting the solar gravity, llnaba and Ikomal (120031 ) derived a re- 
lation between the enhanced radius of a protoplanet and the size of planetesimals, which also 
depends on the atmospheric density distribution (their Eq.(17)). Substituting t he obtained 



enhanced radius into the semi-analytic expressions of collision rates deriv ed bv llnaba et al. 
(I2OOII) . accretion rates for protoplanets with atmospheres are obtained (llnaba and Ikoma 
20031 . their Eqs.(20)-(24)). 



Since (Pace) is a quantity that is integrated with r espect to e, i, an d fe, (P arn) is not 
an explicit function of e, 2, and h. However, (Pace) of llnaba and Ikomal (120031 ) ((Peoi) in 
their notation) is a function of Re-, which explicitly depends on e, z, and h. Here we assume 

that = \j {e*Y -\- for Pc in calculating their (Paee)- In the following comparison, we 
assume an atmospheric structure described by a power-law function (i.e., Eq.(j4|)). 

Figure [TOl shows normalized ac cretion rates obtained by our semi-analytic expression 



and that of 



Inaba and Ikomal (120031 ) when ^ = 10 ^ and r 



lo- 



in this case, all accret- 



ing objects are practically captured by the gas drag of the atmosphere, rather than direct 
collision. We can see that the overall features of the two results are similar to each other; 
the accretion rates in the low velocity regime roughly agree with each other and those in the 
high velocity regime have the same slope. However (Paee) for ^ = 10~^ of our result is about 
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a factor of 2.3 higher than theirs. We find that the difference of the factor of 2.3 between the 
two formulas comes from the following three effects: (i) The amount of energy dissipation 
for in coming particles in an atmosphere is estimated by their Eq.(16) in llnaba and Ikoma 
( 120031 ). This estimate is smaller by a factor of about two compared to the more accurate esti- 
mate in the present work. The correction factor is denoted by fc in the present work and was 
estimated empirically to be 1.83. This factor appears in the expression of the accretion rate 



as jP^r r) oc /c^" (Eq. (l38!l ). thus it enhances (Pace) by ~ 1.83^/"^ = 1.22. (ii) Inaba and Ikoma 
(120031 ) does not seem to fully take account of the 6-dependence of Rc in obtaining the ac- 
cretion rate averaged over the Rayleigh distribution. If it is correctly taken into account, 
this enhances the accretion rate by 1.17. (iii) Instead of fully taking account of the e- and 



i-dependence of Rc and the averaging over the Rayleigh distribution, v^o = y (e*)^ + (z*)^ 
seems to have been used in their calculation of (Pcoi)- If we fully take account of this, the 
accretion rate is enhanc ed by about 1.6. Altho ugh all of their calculation procedures are 
not explicitly written in llnaba and Ikomal (120031 ) . the product of the above three factors is 
about 2.3, which is consistent with the difference observed in Fig. [TD], thus we believe that 
the difference comes from the three factors described above. Based on the consideration, 
we can deduce the reason why the two accretion rates in the low velocity regime roughly 
agree with each other. As described above, two causes out of the three are related to the 
dependence of Rc on e, i and b, while Rc does not depend on them in the low velocity regime. 
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Fig. 1. — Enhanced radius Rc as a function of Voo when a = 3 and = 2. Four hues 
correspond to the cases with ^ = 10~^, 10~^, 10"'', 10~^ from top to bottom. Core radius in 
this normahzation is typically about 0.005 for lAU and 0.001 for 5AU. 
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Fig. 2. — Example orbits to show the effect of gas drag when b = 4.855, e = 3.0, i = 0, 
T = 0.046. Sohd hues show the orbits with gas drag when ^ = 1.66 x 10"^, 1.67 x 10~^, 
2.00 X 10~^, respectively from left to right. Dashed lines in the left and middle panels show 
the orbit without gas drag. Dotted line represents the Hill sphere. 
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Fig. 3. — Change of the Jacobi energy of a particle due to gas drag, as a function of time 
(left panel) and the distance from the origin (right panel), for the orbit shown in the middle 
panel of Figj2] (^ = 1.67 x lO'^). 




Fig. 4. — Total change of the energy E due to gas drag through an encounter with the planet 
for each orbit, as a function of the minimum distance to the planet center in the case with 
a = 3.113. Red pluses and green crosses show the low- velocity (e = 3) and high- velocity 
(e = 30) cases, respectively. For each value of e, numerical results for ,^ = 1 x 10~^ with 
a wide range of initial phase angle r and impact parameter b in the two-dimensional case 
{i = 0) are shown. Blue asterisks show the case without gas drag, that is, the error of the 
numerical integration. Long-dashed line shows the approximate analytic results (Eq.( fT2i) ) 
for the high- velocity regime (i.e., the term proportional to fl~^^) and short-dashed line shows 
that for the low- velocity regime (the term proportional to r~°:J. We set fc = 2. 
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Fig. 5. — Total change of the energy E due to gas drag through an encounter with the 
planet for each orbit, as a function of the minimum distance in the case with a = 3.113. 
Left panel shows the dependence on ^ when {e,i) = (0.1,0.05): red pluses, green crosses, 
and blue asterisks are the cases with ^ = 2.56 x 10~^ (weak), 2.56 x 10~^ (intermediate), 
and 2.56 x 10~^ (strong), respectively. Right panel shows the dependence on initial orbital 
elements when ^ = 2.56 x 10"'': red pluses, green crosses, and blue asterisks are the case 
with {e,i) = (10,5), (1,0.5), and (0.1,0.05), respectively. The vertical dashed line shows the 
planet's radius fp = 5 x 10"'^. 
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Fig. 6. — Differential non-dimensional accretion rate pcap and collision rate pcoi as a function 
of b. Top and bottom panels shows (e, = (17.8,8.9) and (3.16,1.58), respectively. Solid 
and dotted lines show analytic solutions of Pcap and Pcoi, respectively. Pluses and crosses 
show numerical solutions of pcap and pcoi, respectively. 
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Fig. 7. — Normalized accretion rate Pace as a function of e when e = 2i for several values of 
^. Left panel shows the case with fp = 0.005 and right panel shows when fp = 0.001. The 
uppermost line in both panels shows normalized collision probability onto the lemon-shaped 
surface of the Hill sphere. The bottom line in both panels shows Pcoi- 




Fig. 8. — Normalized capture rate as a function of ^ when e = i = 0. Numerical results for 
three values of a are shown by symbols. Six fit lines for these results described by Eqs. fl2Tl) 
and (122|) are also shown. 
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Fig. 9. — Normalized accretion rates averaged over the Rayleigh distribution for fp = 0.005 
(left panel) and fp = 0.001 (right panel), as a function of e*(= 2i*). Marks show simulation 
results for ^ = 10^^, 10^^, 10"^, 10"^, and 10~^, and lines show the semi-analytical formula 
given by Eq. fl26p . The Rayleigh distribution average of collision rate (Eq. (I27|) ) is also shown 
by the long dashed line for comparison in the case of fp = 0.001 (left panel). Note that, 
in the case of Vp = 0.005 and ^ = 10~^, there is no effect of the atmosphere on (Pace), i-e.. 
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Fig. 10. — Normalized accretion rate averaged over the Rayleigh velocity distribution as 
a function of e*(= 2i*) for ^ = 10~^ and fp = 0.001. Solid line shows the accretion rate 



obtained by our ser ni- analytic express i on ('E q. (l26!) ). and dashed hne is obtained by using the 
formulas derived by llnaba and Ikomal (120031 ). In both cases, a power-law density distribution 
is assumed for the atmospheric structure {a = 3). Collision rate for the case without an 
atmosphere is also shown by the dotted line for comparison. 
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Fi g. 11. — Density profile s of a realistic atmosphere derived by the analytic method given 



m 



Inaba and Ikomal (120031 ) for IM^ and lOM^ planets at 5AU (thick-solid lines) in the case 
when pcore = 3.4 x lO^kgm"^, L = 10~^Lq, and = 10~^, and power-law functions fitted 
to the realistic one at the bottom (thick-dashed lines). The exponents for the fit lines are 
a = 3.094 for IM^ and a = 3.027 for IOM0. Thick lines show the case when the temperature 
and gas density at the outer boundary are 125K and 1.5 x 10~^kgm~^, respectively. In order 
to see the dependence on the boundary condition, we also show density profiles when the 
temperature at the boundary is lower (50K, thin dot-dashed line) and the density is lower 
(1.5 X 10"^kgm~^, thin dotted line) in the lOM^ case. 
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Fig. 12. — Change of the energy E due to gas drag as a function of minimum approach 
distance (the same as Fig. [5l) for a reahstic atmosphere model (red pluses) and a power-law 
model (green crosses) when M = lOM^ at 5AU. Left and right panels show the case when 
(e, i) = (0.1,0.05) and (10,5), respectively. Four different planetesimal sizes are assumed for 
each of the two atmosphere models; rg = lO^km, lO^km, 1km, and 10m, from lower-left to 
upper-right. 
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Fig. 13. — Pace for the realistic atmosphere model (thick-solid lines) and for the power-law 
atmosphere model (thick-dotted lines), and Pcoi for = 0.001 (thin dashed line) (e = 2i). 
Left panel shows the cases when M = lOM® with Vg = lO^km (red lines), lO^km (blue lines), 
1km (purple lines), and 10m (green lines). Right panel shows the cases when M = IM® 
with Tg =lkm (purple lines), 10m (green lines), and 0.1m (yellow lines). 
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Fig. 14. — (a) Growth of a protoplanet at 5AU from the Sun in a disk with initial sohd surface 
density of 100kg and with gas density of 3.9 x 10~^kgm~'^. Planetesimals' material 
density is 2 x lO'^kgm"^. Thin solid line shows the case of 10-km-diameter planetesimals 
and a protoplanet without an atmosphere, an d the thin-dashed l ine shows the case with the 
effect of atmosphere based on the formula of IChambersI (j2006al ) for planetesimals of 10km 
diameter. Thick solid, dashed, and dotted lines show the results based on our formula of the 
accretion rate with the effect of atmosphere, in the case with planetesimal of 10km, 1km, and 
100m, respectively, (b) Time evolution of C, in the three cases with our formula of accretion 
rates for planetesimals with 10km (solid line), 1km (dashed line), and 100m (dotted line), 
respectively. 



